/*
  ivl_funcs.c  -  implementation of mathematical elementary functions

  This file is part of the IVL Library (LIBIVL).
  
  LIBIVL - library for interval computations 
  
  Copyright (c) 2007 Mateusz Paprocki

  The IVL Library is free software; you can redistribute it and/or modify
  it under the terms of the GNU Lesser General Public License as published
  by the Free Software Foundation; either version 2.1 of the License, or
  (at your option) any later version.

  The IVL Library is distributed in the hope that it will be useful, but
  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  License for more details.

  You should have received a copy of the GNU Lesser General Public License
  along with the IVL Library; see the file LICENCE.  If not, write to
  the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
  Boston, MA 02110-1301, USA.
*/

#include "interval.h"

int ivl_signum(ivl_ptr ivl, ivl_srcptr op) {
	int inexact = 0;
	
	if (DOM_UNDEFINED(op))
		ivl->domain = IVL_UNDEFINED;
	else {		
		switch (mpfr_sgn(IVL_LEFT(op))) {		
			case -1:
				switch (mpfr_sgn(IVL_LEFT(op))) {
					case -1:
						inexact = ivl_set_s(ivl, -1);
						ivl->domain = op->domain;
						break;
					case 0:
						inexact = ivl_make_u(ivl, -1, 0);
						ivl->domain = IVL_PARTIAL;
						break;
					case 1:
						inexact = ivl_make_u(ivl, -1, 1);
						ivl->domain = IVL_PARTIAL;
						break;
				}
				
				break;
			case 0:
				if (mpfr_sgn(IVL_RIGHT(op)) == 0) {								
					inexact = ivl_set_u(ivl, 1);
					ivl->domain = op->domain;
				} else {
					inexact = ivl_make_u(ivl, 0, 1);
					ivl->domain = IVL_PARTIAL; 
				}
				
				break;
			case 1:
				inexact = ivl_set_u(ivl, 1);
				ivl->domain = op->domain;
				break;
		}
	}
	
	return inexact;
}

int ivl_neg(ivl_ptr ivl, ivl_srcptr op) {
	int inexact = 0, inexact_left, inexact_right;
	
	if (DOM_UNDEFINED(op))
		ivl->domain = IVL_UNDEFINED;
	else {
		/* negation is floating-point exact, but possibly precision is different */ 
		inexact_left = mpfr_neg(IVL_LEFT(ivl), IVL_RIGHT(ivl), GMP_RNDD);
		inexact_right = mpfr_neg(IVL_RIGHT(ivl), IVL_LEFT(ivl), GMP_RNDU);
		
		IVL_CHECK(ivl, op->domain, inexact, inexact_left, inexact_right);						
	}
	
	return inexact;
}

int ivl_abs(ivl_ptr ivl, ivl_srcptr op) {
	int inexact = 0, inexact_left, inexact_right;
	
	if (DOM_UNDEFINED(op))
		ivl->domain = IVL_UNDEFINED;
	else {
		if (mpfr_sgn(IVL_RIGHT(op)) >= 0) {					
			if (mpfr_sgn(IVL_LEFT(op)) >= 0) {
				inexact_left = mpfr_abs(IVL_LEFT(ivl), IVL_LEFT(op), GMP_RNDD);
				inexact_right = mpfr_abs(IVL_RIGHT(ivl), IVL_RIGHT(op), GMP_RNDU);
			} else {
				inexact_left = mpfr_set_ui(IVL_LEFT(ivl), 0, GMP_RNDD);
				
				if (mpfr_cmpabs(IVL_LEFT(op), IVL_RIGHT(op)) > 0)
					inexact_right = mpfr_abs(IVL_RIGHT(ivl), IVL_LEFT(op), GMP_RNDU);
				else
					inexact_right = mpfr_abs(IVL_RIGHT(ivl), IVL_RIGHT(op), GMP_RNDU);
			}			
		} else {
			inexact_left = mpfr_abs(IVL_LEFT(ivl), IVL_RIGHT(op), GMP_RNDD);
			inexact_right = mpfr_abs(IVL_RIGHT(ivl), IVL_LEFT(op), GMP_RNDU);
		}
		
		IVL_CHECK(ivl, op->domain, inexact, inexact_left, inexact_right);		
	}
	
	return inexact;
}

int ivl_floor(ivl_ptr ivl, ivl_srcptr op) {
	int inexact = 0, inexact_left, inexact_right;
	
	if (DOM_UNDEFINED(op))
		ivl->domain = IVL_UNDEFINED;
	else {		
		inexact_left = mpfr_floor(IVL_LEFT(ivl), IVL_LEFT(op));
		inexact_right = mpfr_floor(IVL_RIGHT(ivl), IVL_RIGHT(op));
		
		if (DOM_PARTIAL(op))
			ivl->domain = IVL_PARTIAL;
		else {
			if (mpfr_cmp(IVL_LEFT(ivl), IVL_RIGHT(ivl)))
				ivl->domain = IVL_PARTIAL;
			else
				ivl->domain = IVL_TOTAL;
		}
		
		IVL_CHECK_NO_DOM(ivl, inexact, inexact_left, inexact_right);
	}
	
	return inexact;
}

int ivl_ceil(ivl_ptr ivl, ivl_srcptr op) {
	int inexact = 0, inexact_left, inexact_right;
	
	if (DOM_UNDEFINED(op))
		ivl->domain = IVL_UNDEFINED;
	else {
		inexact_left = mpfr_ceil(IVL_LEFT(ivl), IVL_LEFT(op));
		inexact_right = mpfr_ceil(IVL_RIGHT(ivl), IVL_RIGHT(op));
		
		if (DOM_PARTIAL(op))
			ivl->domain = IVL_PARTIAL;
		else {
			if (mpfr_cmp(IVL_LEFT(ivl), IVL_RIGHT(ivl)))
				ivl->domain = IVL_PARTIAL;
			else
				ivl->domain = IVL_TOTAL;
		}
		
		IVL_CHECK_NO_DOM(ivl, inexact, inexact_left, inexact_right);
	}
	
	return inexact;	
}

int ivl_sqrt(ivl_ptr ivl, ivl_srcptr op) {
	int inexact = 0, inexact_left, inexact_right;
	
	if (DOM_UNDEFINED(op) || (mpfr_sgn(IVL_RIGHT(op)) < 0))
		ivl->domain = IVL_UNDEFINED;
	else {				
		if (mpfr_sgn(IVL_LEFT(op)) >= 0) {
			inexact_left = mpfr_sqrt(IVL_LEFT(ivl), IVL_LEFT(op), GMP_RNDD);	
			ivl->domain = op->domain;
		} else {
			inexact_left = mpfr_set_ui(IVL_LEFT(ivl), 0, GMP_RNDD);	
			ivl->domain = IVL_PARTIAL;
		}
		
		inexact_right = mpfr_sqrt(IVL_RIGHT(ivl), IVL_RIGHT(op), GMP_RNDU);

		IVL_CHECK_NO_DOM(ivl, inexact, inexact_left, inexact_right);
	}
	
	return inexact;
}

int ivl_exp(ivl_ptr ivl, ivl_srcptr op) {
	int inexact = 0, inexact_left, inexact_right;
	
	if (DOM_UNDEFINED(op))
		ivl->domain = IVL_UNDEFINED;
	else {	
		inexact_left = mpfr_exp(IVL_LEFT(ivl), IVL_LEFT(op), GMP_RNDD);
		inexact_right = mpfr_exp(IVL_RIGHT(ivl), IVL_RIGHT(op), GMP_RNDU);

		IVL_CHECK(ivl, op->domain, inexact, inexact_left, inexact_right);
	}
	
	return inexact;	
}

int ivl_exp2(ivl_ptr ivl, ivl_srcptr op) {
	int inexact = 0, inexact_left, inexact_right;
	
	if (DOM_UNDEFINED(op))
		ivl->domain = IVL_UNDEFINED;
	else {	
		inexact_left = mpfr_exp2(IVL_LEFT(ivl), IVL_LEFT(op), GMP_RNDD);
		inexact_right = mpfr_exp2(IVL_RIGHT(ivl), IVL_RIGHT(op), GMP_RNDU);

		IVL_CHECK(ivl, op->domain, inexact, inexact_left, inexact_right);
	}
	
	return inexact;	
}

int ivl_exp10(ivl_ptr ivl, ivl_srcptr op) {
	int inexact = 0, inexact_left, inexact_right;
	
	if (DOM_UNDEFINED(op))
		ivl->domain = IVL_UNDEFINED;
	else {	
		inexact_left = mpfr_exp10(IVL_LEFT(ivl), IVL_LEFT(op), GMP_RNDD);
		inexact_right = mpfr_exp10(IVL_RIGHT(ivl), IVL_RIGHT(op), GMP_RNDU);

		IVL_CHECK(ivl, op->domain, inexact, inexact_left, inexact_right);
	}
	
	return inexact;	
}

int ivl_log(ivl_ptr ivl, ivl_srcptr op) {
	int inexact = 0, inexact_left, inexact_right;
	
	if (DOM_UNDEFINED(op))
		ivl->domain = IVL_UNDEFINED;
	else {
		if (mpfr_sgn(IVL_RIGHT(op)) < 0)
			ivl->domain = IVL_UNDEFINED;			
		else {
			inexact_right = mpfr_log(IVL_RIGHT(ivl), IVL_RIGHT(op), GMP_RNDU);
			
			if (mpfr_sgn(IVL_LEFT(op)) >= 0) {
				inexact_left = mpfr_log(IVL_LEFT(ivl), IVL_LEFT(op), GMP_RNDD);
				ivl->domain = op->domain;
			} else {
				mpfr_set_inf(IVL_LEFT(ivl), -1);
				ivl->domain = IVL_PARTIAL;
				inexact_left = 0;			
			}
	
			IVL_CHECK_NO_DOM(ivl, inexact, inexact_left, inexact_right);
		}
	}
	
	return inexact;
}

int ivl_log2(ivl_ptr ivl, ivl_srcptr op) {
	int inexact = 0, inexact_left, inexact_right;
	
	if (DOM_UNDEFINED(op))
		ivl->domain = IVL_UNDEFINED;
	else {
		if (mpfr_sgn(IVL_RIGHT(op)) < 0)
			ivl->domain = IVL_UNDEFINED;			
		else {
			inexact_right = mpfr_log2(IVL_RIGHT(ivl), IVL_RIGHT(op), GMP_RNDU);
			
			if (mpfr_sgn(IVL_LEFT(op)) >= 0) {
				inexact_left = mpfr_log2(IVL_LEFT(ivl), IVL_LEFT(op), GMP_RNDD);
				ivl->domain = op->domain;
			} else {
				mpfr_set_inf(IVL_LEFT(ivl), -1);
				ivl->domain = IVL_PARTIAL;
				inexact_left = 0;			
			}
	
			IVL_CHECK_NO_DOM(ivl, inexact, inexact_left, inexact_right);
		}
	}
	
	return inexact;
}

int ivl_log10(ivl_ptr ivl, ivl_srcptr op) {
	int inexact = 0, inexact_left, inexact_right;
	
	if (DOM_UNDEFINED(op))
		ivl->domain = IVL_UNDEFINED;
	else {
		if (mpfr_sgn(IVL_RIGHT(op)) < 0)
			ivl->domain = IVL_UNDEFINED;			
		else {
			inexact_right = mpfr_log10(IVL_RIGHT(ivl), IVL_RIGHT(op), GMP_RNDU);
			
			if (mpfr_sgn(IVL_LEFT(op)) >= 0) {
				inexact_left = mpfr_log10(IVL_LEFT(ivl), IVL_LEFT(op), GMP_RNDD);
				ivl->domain = op->domain;
			} else {
				mpfr_set_inf(IVL_LEFT(ivl), -1);
				ivl->domain = IVL_PARTIAL;
				inexact_left = 0;			
			}
	
			IVL_CHECK_NO_DOM(ivl, inexact, inexact_left, inexact_right);
		}
	}
	
	return inexact;
}

int ivl_min(ivl_ptr ivl, ivl_srcptr x, ivl_srcptr y) {
	int inexact = 0, inexact_left, inexact_right;
	
	if (DOM_UNDEFINED(x) || DOM_UNDEFINED(y))
		ivl->domain = IVL_UNDEFINED;
	else {
		if (mpfr_lessequal_p(IVL_LEFT(x), IVL_LEFT(y))) {
			if (mpfr_greater_p(IVL_RIGHT(x), IVL_LEFT(y))) {
				inexact_right = mpfr_set(IVL_RIGHT(ivl), IVL_LEFT(y), GMP_RNDU);
				ivl->domain = x->domain | y->domain;
			} else
				ivl->domain = x->domain;
				
			inexact_left = mpfr_set(IVL_LEFT(ivl), IVL_LEFT(x), GMP_RNDD);
		} else {
			if (mpfr_greater_p(IVL_RIGHT(y), IVL_LEFT(x))) {
				inexact_right = mpfr_set(IVL_RIGHT(ivl), IVL_LEFT(x), GMP_RNDU);
				ivl->domain = x->domain | y->domain;
			} else
				ivl->domain = y->domain;
			
			inexact_left = mpfr_set(IVL_LEFT(ivl), IVL_LEFT(y), GMP_RNDD);
		}
		
		IVL_CHECK_NO_DOM(ivl, inexact, inexact_left, inexact_right);
	}
	
	return inexact;
}

int ivl_max(ivl_ptr ivl, ivl_srcptr x, ivl_srcptr y) {
	int inexact = 0, inexact_left, inexact_right;
	
	if (DOM_UNDEFINED(x) || DOM_UNDEFINED(y))
		ivl->domain = IVL_UNDEFINED;
	else {
		if (mpfr_greaterequal_p(IVL_RIGHT(x), IVL_RIGHT(y))) {
			if (mpfr_less_p(IVL_LEFT(x), IVL_RIGHT(y))) {
				inexact_left = mpfr_set(IVL_LEFT(ivl), IVL_RIGHT(y), GMP_RNDD);
				ivl->domain = x->domain | y->domain;
			} else
				ivl->domain = x->domain;
				
			inexact_right = mpfr_set(IVL_RIGHT(ivl), IVL_RIGHT(x), GMP_RNDU);
		} else {
			if (mpfr_less_p(IVL_LEFT(y), IVL_RIGHT(x))) {
				inexact_left = mpfr_set(IVL_LEFT(ivl), IVL_RIGHT(y), GMP_RNDD);
				ivl->domain = x->domain | y->domain;
			} else
				ivl->domain = y->domain;
			
			inexact_right = mpfr_set(IVL_RIGHT(ivl), IVL_RIGHT(y), GMP_RNDU);
		}
		
		IVL_CHECK_NO_DOM(ivl, inexact, inexact_left, inexact_right);
	}
	
	return inexact;
}

